### Copy thy neighbour: Spatial interdependences in the democracy-repression nexus ###
### Journal of Human Rights ###
### Created by Roman-Gabriel Olar on 01/05/2019###
### Last modified on 25/08/2023###

### Load packages
require(foreign)
require(sandwich)
require(lmtest)
require(MASS)
require(Hmisc)
require(stargazer)
require(sciplot)
require(lme4)

### Clear space
rm(list = ls())


### Define function to combine multiple model estimates (this incorporates uncertainty from the latent variables)
milm <- function(fml, midata){
  xx <- terms(as.formula(fml))
  lms <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
  ses <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
  vcovs <- list()
  for(i in 1:length(midata)){
    tmp <- lm(formula=as.formula(fml), data=midata[[i]])
    lms[,i] <- tmp$coefficients
    ses[,i] <- sqrt(diag(vcov(tmp)))
    vcovs[[i]] <- vcov(tmp)
  }
  par.est <- apply(lms, 1, mean)
  se.within <- apply(ses, 1, mean)
  se.between <- apply(lms, 1, var)
  se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(midata))))
  list("terms"=names(tmp$coefficients), "beta" = par.est, "SE"=se.est, "vcovs"=vcovs,"coefs" = lms)    
}

### Load data. 
setwd("set your own working directory")

dat <- read.dta("MainData.dta")

attach(dat)
datfh <- na.omit(data.frame(theta_mean, theta_sd, l_fariss, l_fariss_sd, SLreg_fariss_t1, SLreg_fariss_sd, v2cltort_sd, v2clkill_sd,
                            l_polity, l_polyarchy, l_polyarchy_sd, SLreg_polyarchy_t1, SLreg_polyarchy_sd, l_uds, l_boix, l_cgv, l_gwf_democracy, 
                            SLreg_polity, SLreg_gwf, SLreg_polyarchy_t1, SLreg_boix_t1, SLreg_cgv_t1,
                            SLreg_clkill_sd, SLreg_clkill_t1, SLreg_cltort_sd, SLreg_cltort_t1, l_v2clkill_sd, l_v2clkill, l_cltort_sd, l_cltort, v2cltort, v2clkill,
                            l_uds2, l_uds_sd, SLreg_uds_t1, SLreg_uds_sd, l_gdp_gro, l_lngdppc, l_lnpop, civil_war, international_war, l_resdep,
                            ccode_fe1, ccode_fe2, ccode_fe3, ccode_fe4, ccode_fe5, ccode_fe6, ccode_fe7, ccode_fe8, ccode_fe9, 
                            ccode_fe10, ccode_fe11, ccode_fe12, ccode_fe13, ccode_fe14, ccode_fe15, ccode_fe16, ccode_fe17, 
                            ccode_fe18, ccode_fe19, ccode_fe20, ccode_fe21, ccode_fe22, ccode_fe23, ccode_fe24, ccode_fe25, 
                            ccode_fe26, ccode_fe27, ccode_fe28, ccode_fe29, ccode_fe30, ccode_fe31, ccode_fe32, ccode_fe33, 
                            ccode_fe34, ccode_fe35, ccode_fe36, ccode_fe37, ccode_fe38, ccode_fe39, ccode_fe40, ccode_fe41, 
                            ccode_fe42, ccode_fe43, ccode_fe44, ccode_fe45, ccode_fe46, ccode_fe47, ccode_fe48, ccode_fe49, 
                            ccode_fe50, ccode_fe51, ccode_fe52, ccode_fe53, ccode_fe54, ccode_fe55, ccode_fe56, ccode_fe57,
                            ccode_fe58, ccode_fe59, ccode_fe60, ccode_fe61, ccode_fe62, ccode_fe63, ccode_fe64, ccode_fe65, 
                            ccode_fe66, ccode_fe67, ccode_fe68, ccode_fe69, ccode_fe70, ccode_fe71, ccode_fe72, ccode_fe73, 
                            ccode_fe74, ccode_fe75, ccode_fe76, ccode_fe77, ccode_fe78, ccode_fe79, ccode_fe80, ccode_fe81, 
                            ccode_fe82, ccode_fe83, ccode_fe84, ccode_fe85, ccode_fe86, ccode_fe87, ccode_fe88, ccode_fe89, 
                            ccode_fe90, ccode_fe91, ccode_fe92, ccode_fe93, ccode_fe94, ccode_fe95, ccode_fe96, ccode_fe97, 
                            ccode_fe98, ccode_fe99, ccode_fe100, ccode_fe101, ccode_fe102, ccode_fe103, ccode_fe104, ccode_fe105, 
                            ccode_fe106, ccode_fe107, ccode_fe108, ccode_fe109, ccode_fe110, ccode_fe111, ccode_fe112, ccode_fe113, 
                            ccode_fe114, ccode_fe115, ccode_fe116, ccode_fe117, ccode_fe118, ccode_fe119, ccode_fe120, ccode_fe121, 
                            ccode_fe122, ccode_fe123, ccode_fe124, ccode_fe125, ccode_fe126, ccode_fe127, ccode_fe128, ccode_fe129, 
                            ccode_fe130, ccode_fe131, ccode_fe132, ccode_fe133, ccode_fe134, ccode_fe135, ccode_fe136, ccode_fe137, 
                            ccode_fe138, y_fe1, y_fe2, y_fe3, y_fe4, y_fe5, y_fe6, 
                            y_fe7, y_fe8, y_fe9, y_fe10, y_fe11, y_fe12, y_fe13, y_fe14, y_fe15, y_fe16, y_fe17, y_fe18, y_fe19, y_fe20, 
                            y_fe21, y_fe22, y_fe23, y_fe24, y_fe25, y_fe26, y_fe27, y_fe28, y_fe29, y_fe30, y_fe31, y_fe32, y_fe33, y_fe34, 
                            y_fe35, y_fe36, y_fe37, y_fe38, y_fe39, y_fe40, y_fe41, y_fe42, y_fe43, y_fe44, y_fe45, y_fe46, y_fe47, y_fe48, 
                            y_fe49, y_fe50, y_fe51, y_fe52, y_fe53, y_fe54, y_fe55, y_fe56, y_fe57, y_fe58, y_fe59, y_fe60, y_fe61))

dim(datfh)


### Setup data to run models where we account for uncertainty
set.seed(32261991)
newdata <- list()
for(i in 1:1000){
  newdata[[i]] <- datfh
  newdata[[i]]$draw <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$theta_mean, sd=datfh$theta_sd)
  newdata[[i]]$draw.lag <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$l_fariss, sd=datfh$l_fariss_sd)
  newdata[[i]]$draw.lag.SLfariss <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$SLreg_fariss_t1, sd=datfh$SLreg_fariss_sd)
  newdata[[i]]$draw.lag.SLpolyarchy <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$SLreg_polyarchy_t1, sd=datfh$SLreg_polyarchy_sd)
  newdata[[i]]$draw.lag.polyarchy <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$l_polyarchy, sd=datfh$l_polyarchy_sd)
}

### Table 3 Analysis ###

### Model 1 - Non-spatial (NS)
FML_NS_VDEM <-"draw ~ draw.lag + draw.lag.polyarchy +
l_gdp_gro + l_lngdppc + l_lnpop + civil_war + international_war + l_resdep"

model_ns_vdem <- milm(fml=FML_NS_VDEM, midata=newdata)
df <- 5870

# create object for printing during the loop and then print the results to screen
model1<- data.frame(model_ns_vdem$terms, model_ns_vdem$beta, model_ns_vdem$SE, model_ns_vdem$beta/model_ns_vdem$SE)
colnames(model1) <- c("Variable","Est. Coef.","SE","t-statistic")

print(paste("Model Output"))
print(model1)


### Model 2 - Non-Spatial + FE
FML_NS_FE_VDEM <-"draw ~ draw.lag + draw.lag.polyarchy +
l_gdp_gro + l_lngdppc + l_lnpop + civil_war + international_war + l_resdep + ccode_fe1+ 
ccode_fe2+ ccode_fe3+ ccode_fe4 + ccode_fe5+ ccode_fe6+ ccode_fe7+ ccode_fe8+ ccode_fe9+ 
ccode_fe10+ ccode_fe11+ ccode_fe12+ ccode_fe13+ ccode_fe14+ ccode_fe15+ ccode_fe16+ ccode_fe17+ 
ccode_fe18+ ccode_fe19+ ccode_fe20+ ccode_fe21+ ccode_fe22+ ccode_fe23+ ccode_fe24+ ccode_fe25+ 
ccode_fe26+ ccode_fe27+ ccode_fe28+ ccode_fe29+ ccode_fe30+ ccode_fe31+ ccode_fe32+ ccode_fe33+ 
ccode_fe34+ ccode_fe35+ ccode_fe36+ ccode_fe37+ ccode_fe38+ ccode_fe39+ ccode_fe40+ ccode_fe41+ 
ccode_fe42+ ccode_fe43+ ccode_fe44+ ccode_fe45+ ccode_fe46+ ccode_fe47+ ccode_fe48+ ccode_fe49+ 
ccode_fe50+ ccode_fe51+ ccode_fe52+ ccode_fe53+ ccode_fe54+ ccode_fe55+ ccode_fe56+ ccode_fe57+
ccode_fe58+ ccode_fe59+ ccode_fe60+ ccode_fe61+ ccode_fe62+ ccode_fe63+ ccode_fe64+ ccode_fe65+ 
ccode_fe66+ ccode_fe67+ ccode_fe68+ ccode_fe69+ ccode_fe70+ ccode_fe71+ ccode_fe72+ ccode_fe73+ 
ccode_fe74+ ccode_fe75+ ccode_fe76+ ccode_fe77+ ccode_fe78+ ccode_fe79+ ccode_fe80+ ccode_fe81+ 
ccode_fe82+ ccode_fe83+ ccode_fe85+ ccode_fe86+ ccode_fe87+ ccode_fe88+ ccode_fe89+ 
ccode_fe90+ ccode_fe91+ ccode_fe92+ ccode_fe93+ ccode_fe94+ ccode_fe95+ ccode_fe96+ ccode_fe97+ 
ccode_fe98+ ccode_fe99+ ccode_fe100+ ccode_fe101+ ccode_fe102+ ccode_fe103+ ccode_fe104+ ccode_fe105+ 
ccode_fe106+ ccode_fe107+ ccode_fe108+ ccode_fe109+ ccode_fe110+ ccode_fe111+ ccode_fe112+ ccode_fe113+ 
ccode_fe114+ ccode_fe115+ ccode_fe116+ ccode_fe117+ ccode_fe118+ ccode_fe119+ ccode_fe120+ ccode_fe121+ 
ccode_fe122+ ccode_fe123+ ccode_fe124+ ccode_fe125+ ccode_fe126+ ccode_fe127+ ccode_fe128+ ccode_fe129+ 
ccode_fe130+ ccode_fe131+ ccode_fe132+ ccode_fe133+ ccode_fe134+ ccode_fe135+ ccode_fe136+ ccode_fe137+ 
ccode_fe138+ y_fe1+ y_fe2+ y_fe3+ y_fe5+ y_fe6+ 
y_fe7+ y_fe8+ y_fe9+ y_fe10+ y_fe11+ y_fe12+ y_fe13+ y_fe14+ y_fe15+ y_fe16+ y_fe17+ y_fe18+ y_fe19+ y_fe20+ 
y_fe21+ y_fe22+ y_fe23+ y_fe24+ y_fe25+ y_fe26+ y_fe27+ y_fe28+ y_fe29+ y_fe30+ y_fe31+ y_fe32+ y_fe33+ y_fe34+ 
y_fe35+ y_fe36+ y_fe37+ y_fe38+ y_fe39+ y_fe40+ y_fe41+ y_fe42+ y_fe43+ y_fe44+ y_fe45+ y_fe46+ y_fe47+ y_fe48+ 
y_fe49+ y_fe50+ y_fe51+ y_fe52+ y_fe53+ y_fe54+ y_fe55+ y_fe56+ y_fe57+ y_fe58+ y_fe59+ y_fe60+ y_fe61"


model_ns_fe_vdem <- milm(fml=FML_NS_FE_VDEM, midata=newdata)

model2<- data.frame(model_ns_fe_vdem$terms, model_ns_fe_vdem$beta, model_ns_fe_vdem$SE, model_ns_fe_vdem$beta/model_ns_fe_vdem$SE)
colnames(model2) <- c("Variable","Est. Coef.","SE","z-score")
print(paste("V-Dem NS + FE"))
print(model2)

### Model 3 - Spatial Autoregressive (SAR)
FML_SAR_VDEM <-"draw ~ draw.lag + draw.lag.SLfariss + draw.lag.polyarchy +
l_gdp_gro + l_lngdppc + l_lnpop + civil_war + international_war + l_resdep + ccode_fe1+ 
ccode_fe2+ ccode_fe3+ ccode_fe4 + ccode_fe5+ ccode_fe6+ ccode_fe7+ ccode_fe8+ ccode_fe9+ 
ccode_fe10+ ccode_fe11+ ccode_fe12+ ccode_fe13+ ccode_fe14+ ccode_fe15+ ccode_fe16+ ccode_fe17+ 
ccode_fe18+ ccode_fe19+ ccode_fe20+ ccode_fe21+ ccode_fe22+ ccode_fe23+ ccode_fe24+ ccode_fe25+ 
ccode_fe26+ ccode_fe27+ ccode_fe28+ ccode_fe29+ ccode_fe30+ ccode_fe31+ ccode_fe32+ ccode_fe33+ 
ccode_fe34+ ccode_fe35+ ccode_fe36+ ccode_fe37+ ccode_fe38+ ccode_fe39+ ccode_fe40+ ccode_fe41+ 
ccode_fe42+ ccode_fe43+ ccode_fe44+ ccode_fe45+ ccode_fe46+ ccode_fe47+ ccode_fe48+ ccode_fe49+ 
ccode_fe50+ ccode_fe51+ ccode_fe52+ ccode_fe53+ ccode_fe54+ ccode_fe55+ ccode_fe56+ ccode_fe57+
ccode_fe58+ ccode_fe59+ ccode_fe60+ ccode_fe61+ ccode_fe62+ ccode_fe63+ ccode_fe64+ ccode_fe65+ 
ccode_fe66+ ccode_fe67+ ccode_fe68+ ccode_fe69+ ccode_fe70+ ccode_fe71+ ccode_fe72+ ccode_fe73+ 
ccode_fe74+ ccode_fe75+ ccode_fe76+ ccode_fe77+ ccode_fe78+ ccode_fe79+ ccode_fe80+ ccode_fe81+ 
ccode_fe82+ ccode_fe83+ ccode_fe85+ ccode_fe86+ ccode_fe87+ ccode_fe88+ ccode_fe89+ 
ccode_fe90+ ccode_fe91+ ccode_fe92+ ccode_fe93+ ccode_fe94+ ccode_fe95+ ccode_fe96+ ccode_fe97+ 
ccode_fe98+ ccode_fe99+ ccode_fe100+ ccode_fe101+ ccode_fe102+ ccode_fe103+ ccode_fe104+ ccode_fe105+ 
ccode_fe106+ ccode_fe107+ ccode_fe108+ ccode_fe109+ ccode_fe110+ ccode_fe111+ ccode_fe112+ ccode_fe113+ 
ccode_fe114+ ccode_fe115+ ccode_fe116+ ccode_fe117+ ccode_fe118+ ccode_fe119+ ccode_fe120+ ccode_fe121+ 
ccode_fe122+ ccode_fe123+ ccode_fe124+ ccode_fe125+ ccode_fe126+ ccode_fe127+ ccode_fe128+ ccode_fe129+ 
ccode_fe130+ ccode_fe131+ ccode_fe132+ ccode_fe133+ ccode_fe134+ ccode_fe135+ ccode_fe136+ ccode_fe137+ 
ccode_fe138+ y_fe1+ y_fe2+ y_fe3+ y_fe5+ y_fe6+ 
y_fe7+ y_fe8+ y_fe9+ y_fe10+ y_fe11+ y_fe12+ y_fe13+ y_fe14+ y_fe15+ y_fe16+ y_fe17+ y_fe18+ y_fe19+ y_fe20+ 
y_fe21+ y_fe22+ y_fe23+ y_fe24+ y_fe25+ y_fe26+ y_fe27+ y_fe28+ y_fe29+ y_fe30+ y_fe31+ y_fe32+ y_fe33+ y_fe34+ 
y_fe35+ y_fe36+ y_fe37+ y_fe38+ y_fe39+ y_fe40+ y_fe41+ y_fe42+ y_fe43+ y_fe44+ y_fe45+ y_fe46+ y_fe47+ y_fe48+ 
y_fe49+ y_fe50+ y_fe51+ y_fe52+ y_fe53+ y_fe54+ y_fe55+ y_fe56+ y_fe57+ y_fe58+ y_fe59+ y_fe60+ y_fe61"

model_sar_vdem <- milm(fml=FML_SAR_VDEM, midata=newdata)

model3<- data.frame(model_sar_vdem$terms, model_sar_vdem$beta, model_sar_vdem$SE, model_sar_vdem$beta/model_sar_vdem$SE)
colnames(model3) <- c("Variable","Est. Coef.","SE","z-score")
print(paste("SAR V-Dem"))
print(model3)


### Model 4 - Spatial Lag of X (SLX)

FML_SLX_VDEM <-"draw ~ draw.lag + draw.lag.SLpolyarchy + draw.lag.polyarchy +
l_gdp_gro + l_lngdppc + l_lnpop + civil_war + international_war + l_resdep + ccode_fe1+ 
ccode_fe2+ ccode_fe3+ ccode_fe4 + ccode_fe5+ ccode_fe6+ ccode_fe7+ ccode_fe8+ ccode_fe9+ 
ccode_fe10+ ccode_fe11+ ccode_fe12+ ccode_fe13+ ccode_fe14+ ccode_fe15+ ccode_fe16+ ccode_fe17+ 
ccode_fe18+ ccode_fe19+ ccode_fe20+ ccode_fe21+ ccode_fe22+ ccode_fe23+ ccode_fe24+ ccode_fe25+ 
ccode_fe26+ ccode_fe27+ ccode_fe28+ ccode_fe29+ ccode_fe30+ ccode_fe31+ ccode_fe32+ ccode_fe33+ 
ccode_fe34+ ccode_fe35+ ccode_fe36+ ccode_fe37+ ccode_fe38+ ccode_fe39+ ccode_fe40+ ccode_fe41+ 
ccode_fe42+ ccode_fe43+ ccode_fe44+ ccode_fe45+ ccode_fe46+ ccode_fe47+ ccode_fe48+ ccode_fe49+ 
ccode_fe50+ ccode_fe51+ ccode_fe52+ ccode_fe53+ ccode_fe54+ ccode_fe55+ ccode_fe56+ ccode_fe57+
ccode_fe58+ ccode_fe59+ ccode_fe60+ ccode_fe61+ ccode_fe62+ ccode_fe63+ ccode_fe64+ ccode_fe65+ 
ccode_fe66+ ccode_fe67+ ccode_fe68+ ccode_fe69+ ccode_fe70+ ccode_fe71+ ccode_fe72+ ccode_fe73+ 
ccode_fe74+ ccode_fe75+ ccode_fe76+ ccode_fe77+ ccode_fe78+ ccode_fe79+ ccode_fe80+ ccode_fe81+ 
ccode_fe82+ ccode_fe83+ ccode_fe85+ ccode_fe86+ ccode_fe87+ ccode_fe88+ ccode_fe89+ 
ccode_fe90+ ccode_fe91+ ccode_fe92+ ccode_fe93+ ccode_fe94+ ccode_fe95+ ccode_fe96+ ccode_fe97+ 
ccode_fe98+ ccode_fe99+ ccode_fe100+ ccode_fe101+ ccode_fe102+ ccode_fe103+ ccode_fe104+ ccode_fe105+ 
ccode_fe106+ ccode_fe107+ ccode_fe108+ ccode_fe109+ ccode_fe110+ ccode_fe111+ ccode_fe112+ ccode_fe113+ 
ccode_fe114+ ccode_fe115+ ccode_fe116+ ccode_fe117+ ccode_fe118+ ccode_fe119+ ccode_fe120+ ccode_fe121+ 
ccode_fe122+ ccode_fe123+ ccode_fe124+ ccode_fe125+ ccode_fe126+ ccode_fe127+ ccode_fe128+ ccode_fe129+ 
ccode_fe130+ ccode_fe131+ ccode_fe132+ ccode_fe133+ ccode_fe134+ ccode_fe135+ ccode_fe136+ ccode_fe137+ 
ccode_fe138+ y_fe1+ y_fe2+ y_fe3+ y_fe5+ y_fe6+ 
y_fe7+ y_fe8+ y_fe9+ y_fe10+ y_fe11+ y_fe12+ y_fe13+ y_fe14+ y_fe15+ y_fe16+ y_fe17+ y_fe18+ y_fe19+ y_fe20+ 
y_fe21+ y_fe22+ y_fe23+ y_fe24+ y_fe25+ y_fe26+ y_fe27+ y_fe28+ y_fe29+ y_fe30+ y_fe31+ y_fe32+ y_fe33+ y_fe34+ 
y_fe35+ y_fe36+ y_fe37+ y_fe38+ y_fe39+ y_fe40+ y_fe41+ y_fe42+ y_fe43+ y_fe44+ y_fe45+ y_fe46+ y_fe47+ y_fe48+ 
y_fe49+ y_fe50+ y_fe51+ y_fe52+ y_fe53+ y_fe54+ y_fe55+ y_fe56+ y_fe57+ y_fe58+ y_fe59+ y_fe60+ y_fe61"


model_slx_vdem <- milm(fml=FML_SLX_VDEM, midata=newdata)

model4<- data.frame(model_slx_vdem$terms, model_slx_vdem$beta, model_slx_vdem$SE, model_slx_vdem$beta/model_slx_vdem$SE)
colnames(model4) <- c("Variable","Est. Coef.","SE","z-score")
print(paste("SLX V-Dem"))
print(model4)


### Model 5 - Spatial Durbin Model (SDM)
FML_SDM_VDEM <-"draw ~ draw.lag + draw.lag.SLfariss + draw.lag.SLpolyarchy + draw.lag.polyarchy +
l_gdp_gro + l_lngdppc + l_lnpop + civil_war + international_war + l_resdep + ccode_fe1+ 
ccode_fe2+ ccode_fe3+ ccode_fe4 + ccode_fe5+ ccode_fe6+ ccode_fe7+ ccode_fe8+ ccode_fe9+ 
ccode_fe10+ ccode_fe11+ ccode_fe12+ ccode_fe13+ ccode_fe14+ ccode_fe15+ ccode_fe16+ ccode_fe17+ 
ccode_fe18+ ccode_fe19+ ccode_fe20+ ccode_fe21+ ccode_fe22+ ccode_fe23+ ccode_fe24+ ccode_fe25+ 
ccode_fe26+ ccode_fe27+ ccode_fe28+ ccode_fe29+ ccode_fe30+ ccode_fe31+ ccode_fe32+ ccode_fe33+ 
ccode_fe34+ ccode_fe35+ ccode_fe36+ ccode_fe37+ ccode_fe38+ ccode_fe39+ ccode_fe40+ ccode_fe41+ 
ccode_fe42+ ccode_fe43+ ccode_fe44+ ccode_fe45+ ccode_fe46+ ccode_fe47+ ccode_fe48+ ccode_fe49+ 
ccode_fe50+ ccode_fe51+ ccode_fe52+ ccode_fe53+ ccode_fe54+ ccode_fe55+ ccode_fe56+ ccode_fe57+
ccode_fe58+ ccode_fe59+ ccode_fe60+ ccode_fe61+ ccode_fe62+ ccode_fe63+ ccode_fe64+ ccode_fe65+ 
ccode_fe66+ ccode_fe67+ ccode_fe68+ ccode_fe69+ ccode_fe70+ ccode_fe71+ ccode_fe72+ ccode_fe73+ 
ccode_fe74+ ccode_fe75+ ccode_fe76+ ccode_fe77+ ccode_fe78+ ccode_fe79+ ccode_fe80+ ccode_fe81+ 
ccode_fe82+ ccode_fe83+ ccode_fe85+ ccode_fe86+ ccode_fe87+ ccode_fe88+ ccode_fe89+ 
ccode_fe90+ ccode_fe91+ ccode_fe92+ ccode_fe93+ ccode_fe94+ ccode_fe95+ ccode_fe96+ ccode_fe97+ 
ccode_fe98+ ccode_fe99+ ccode_fe100+ ccode_fe101+ ccode_fe102+ ccode_fe103+ ccode_fe104+ ccode_fe105+ 
ccode_fe106+ ccode_fe107+ ccode_fe108+ ccode_fe109+ ccode_fe110+ ccode_fe111+ ccode_fe112+ ccode_fe113+ 
ccode_fe114+ ccode_fe115+ ccode_fe116+ ccode_fe117+ ccode_fe118+ ccode_fe119+ ccode_fe120+ ccode_fe121+ 
ccode_fe122+ ccode_fe123+ ccode_fe124+ ccode_fe125+ ccode_fe126+ ccode_fe127+ ccode_fe128+ ccode_fe129+ 
ccode_fe130+ ccode_fe131+ ccode_fe132+ ccode_fe133+ ccode_fe134+ ccode_fe135+ ccode_fe136+ ccode_fe137+ 
ccode_fe138+ y_fe1+ y_fe2+ y_fe3+ y_fe5+ y_fe6+ 
y_fe7+ y_fe8+ y_fe9+ y_fe10+ y_fe11+ y_fe12+ y_fe13+ y_fe14+ y_fe15+ y_fe16+ y_fe17+ y_fe18+ y_fe19+ y_fe20+ 
y_fe21+ y_fe22+ y_fe23+ y_fe24+ y_fe25+ y_fe26+ y_fe27+ y_fe28+ y_fe29+ y_fe30+ y_fe31+ y_fe32+ y_fe33+ y_fe34+ 
y_fe35+ y_fe36+ y_fe37+ y_fe38+ y_fe39+ y_fe40+ y_fe41+ y_fe42+ y_fe43+ y_fe44+ y_fe45+ y_fe46+ y_fe47+ y_fe48+ 
y_fe49+ y_fe50+ y_fe51+ y_fe52+ y_fe53+ y_fe54+ y_fe55+ y_fe56+ y_fe57+ y_fe58+ y_fe59+ y_fe60+ y_fe61"


model_sdm_vdem <- milm(fml=FML_SDM_VDEM, midata=newdata)

model5<- data.frame(model_sdm_vdem$terms, model_sdm_vdem$beta, model_sdm_vdem$SE, model_sdm_vdem$beta/model_sdm_vdem$SE)
colnames(model5) <- c("Variable","Est. Coef.","SE","z-score")
print(paste("SDM V-Dem"))
print(model5)

###PRINT RESULTS###
library(stargazer)
stargazer(model1, type = "html", style = "apsr",
          summary = FALSE, single.row = TRUE, out="m1_vdem.doc")

stargazer(model2, type = "html", style = "apsr",
          summary = FALSE, single.row = TRUE, out="m2_vdem.doc")

stargazer(model3, type = "html", style = "apsr", 
          summary = FALSE, single.row = TRUE, out="m3_vdem.doc")

stargazer(model4, type = "html", style = "apsr", 
          summary = FALSE, single.row = TRUE, out="m4_vdem.doc")

stargazer(model5, type = "html", style = "apsr", 
          summary = FALSE, single.row = TRUE, out="m5_vdem.doc")

### Table 4 Analysis ###

### Clear space
rm(list = ls())

### Load data
dat <- read.dta("mSTAR.dta")

attach(dat)
datfh <- na.omit(data.frame(theta_mean, theta_sd, l_fariss, l_fariss_sd, SLreg_fariss_t1, SLreg_fariss_sd, v2cltort_sd, v2clkill_sd,
                            l_polity, l_polyarchy, l_polyarchy_sd, SLreg_polyarchy_t1, SLreg_polyarchy_sd, l_uds, l_boix, l_cgv, l_gwf_democracy, 
                            SLreg_polity, SLreg_gwf, SLreg_polyarchy_t1, SLreg_boix_t1, SLreg_cgv_t1,
                            SL_alli_reg_fariss, SL_alli_reg_fariss_sd, SL_trade_reg_fariss, SL_trade_reg_fariss_sd, 
                            SL_disc_reg_fariss, SL_disc_reg_fariss_sd, SL_conflict_reg_fariss, SL_conflict_reg_fariss_sd,
                            SL_igo_reg_fariss, SL_igo_reg_fariss_sd, SLcont_fariss, SLcont_fariss_sd,
                            l_uds2, l_uds_sd, SLreg_uds_t1, SLreg_uds_sd, l_gdp_gro, l_lngdppc, l_lnpop, civil_war, international_war, l_resdep,
                            ccode_fe1, ccode_fe2, ccode_fe3, ccode_fe4, ccode_fe5, ccode_fe6, ccode_fe7, ccode_fe8, ccode_fe9, 
                            ccode_fe10, ccode_fe11, ccode_fe12, ccode_fe13, ccode_fe14, ccode_fe15, ccode_fe16, ccode_fe17, 
                            ccode_fe18, ccode_fe19, ccode_fe20, ccode_fe21, ccode_fe22, ccode_fe23, ccode_fe24, ccode_fe25, 
                            ccode_fe26, ccode_fe27, ccode_fe28, ccode_fe29, ccode_fe30, ccode_fe31, ccode_fe32, ccode_fe33, 
                            ccode_fe34, ccode_fe35, ccode_fe36, ccode_fe37, ccode_fe38, ccode_fe39, ccode_fe40, ccode_fe41, 
                            ccode_fe42, ccode_fe43, ccode_fe44, ccode_fe45, ccode_fe46, ccode_fe47, ccode_fe48, ccode_fe49, 
                            ccode_fe50, ccode_fe51, ccode_fe52, ccode_fe53, ccode_fe54, ccode_fe55, ccode_fe56, ccode_fe57,
                            ccode_fe58, ccode_fe59, ccode_fe60, ccode_fe61, ccode_fe62, ccode_fe63, ccode_fe64, ccode_fe65, 
                            ccode_fe66, ccode_fe67, ccode_fe68, ccode_fe69, ccode_fe70, ccode_fe71, ccode_fe72, ccode_fe73, 
                            ccode_fe74, ccode_fe75, ccode_fe76, ccode_fe77, ccode_fe78, ccode_fe79, ccode_fe80, ccode_fe81, 
                            ccode_fe82, ccode_fe83, ccode_fe84, ccode_fe85, ccode_fe86, ccode_fe87, ccode_fe88, ccode_fe89, 
                            ccode_fe90, ccode_fe91, ccode_fe92, ccode_fe93, ccode_fe94, ccode_fe95, ccode_fe96, ccode_fe97, 
                            ccode_fe98, ccode_fe99, ccode_fe100, ccode_fe101, ccode_fe102, ccode_fe103, ccode_fe104, ccode_fe105, 
                            ccode_fe106, ccode_fe107, ccode_fe108, ccode_fe109, ccode_fe110, ccode_fe111, ccode_fe112, ccode_fe113, 
                            ccode_fe114, ccode_fe115, ccode_fe116, ccode_fe117, ccode_fe118, ccode_fe119, ccode_fe120, ccode_fe121, 
                            ccode_fe122, ccode_fe123, ccode_fe124, ccode_fe125, ccode_fe126, ccode_fe127, ccode_fe128, ccode_fe129, 
                            ccode_fe130, ccode_fe131, ccode_fe132, ccode_fe133, ccode_fe134, ccode_fe135, ccode_fe136, ccode_fe137, 
                            ccode_fe138, y_fe1, y_fe2, y_fe3, y_fe4, y_fe5, y_fe6, 
                            y_fe7, y_fe8, y_fe9, y_fe10, y_fe11, y_fe12, y_fe13, y_fe14, y_fe15, y_fe16, y_fe17, y_fe18, y_fe19, y_fe20, 
                            y_fe21, y_fe22, y_fe23, y_fe24, y_fe25, y_fe26, y_fe27, y_fe28, y_fe29, y_fe30, y_fe31, y_fe32, y_fe33, y_fe34, 
                            y_fe35, y_fe36, y_fe37, y_fe38, y_fe39, y_fe40, y_fe41, y_fe42, y_fe43, y_fe44, y_fe45, y_fe46, y_fe47, y_fe48, 
                            y_fe49, y_fe50, y_fe51, y_fe52, y_fe53, y_fe54, y_fe55, y_fe56, y_fe57, y_fe58, y_fe59, y_fe60, y_fe61))

dim(datfh)


### Setup data to run models where we account for uncertainty
set.seed(32261991)
newdata <- list()
for(i in 1:1000){
  newdata[[i]] <- datfh
  newdata[[i]]$draw <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$theta_mean, sd=datfh$theta_sd)
  newdata[[i]]$draw.lag <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$l_fariss, sd=datfh$l_fariss_sd)
  newdata[[i]]$draw.lag.SLtrade <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$SL_trade_reg_fariss, sd=datfh$SL_trade_reg_fariss_sd)
  newdata[[i]]$draw.lag.SLalliance <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$SL_alli_reg_fariss, sd=datfh$SL_alli_reg_fariss_sd)
  newdata[[i]]$draw.lag.SLconflict <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$SL_conflict_reg_fariss, sd=datfh$SL_conflict_reg_fariss_sd)
  newdata[[i]]$draw.lag.SLigo <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$SL_igo_reg_fariss, sd=datfh$SL_igo_reg_fariss_sd)
  newdata[[i]]$draw.lag.SLcont <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$SLcont_fariss, sd=datfh$SLcont_fariss_sd)
  newdata[[i]]$draw.lag.polyarchy <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$l_polyarchy, sd=datfh$l_polyarchy_sd)
}


#Model 5 - multi-Spatial Autoregressive (mSTAR)
mSTAR <-"draw ~ draw.lag + draw.lag.SLigo + draw.lag.SLcont +
                  draw.lag.SLtrade + draw.lag.SLalliance +draw.lag.SLconflict +
                  draw.lag.polyarchy + l_gdp_gro + l_lngdppc + l_lnpop + civil_war + international_war + l_resdep + 
                  ccode_fe1+ ccode_fe2+ ccode_fe3+ ccode_fe4+ ccode_fe5+ ccode_fe6+ ccode_fe7+ ccode_fe8+ ccode_fe9+ 
                  ccode_fe10+ ccode_fe11+ ccode_fe12+ ccode_fe13+ ccode_fe14+ ccode_fe15+ ccode_fe16+ ccode_fe17+ 
                  ccode_fe18+ ccode_fe19+ ccode_fe20+ ccode_fe21+ ccode_fe22+ ccode_fe23+ ccode_fe24+ ccode_fe25+ 
                  ccode_fe26+ ccode_fe27+ ccode_fe28+ ccode_fe30+ ccode_fe31+ ccode_fe32+ ccode_fe33+ 
                  ccode_fe34+ ccode_fe35+ ccode_fe36+ ccode_fe37+ ccode_fe38+ ccode_fe39+ ccode_fe40+ ccode_fe41+ 
                  ccode_fe42+ ccode_fe43+ ccode_fe44+ ccode_fe45+ ccode_fe46+ ccode_fe47+ ccode_fe48+ ccode_fe49+ 
                  ccode_fe50+ ccode_fe51+ ccode_fe52+ ccode_fe53+ ccode_fe54+ ccode_fe55+ ccode_fe56+ ccode_fe57+
                  ccode_fe58+ ccode_fe59+ ccode_fe60+ ccode_fe61+ ccode_fe62+ ccode_fe63+ ccode_fe64+ ccode_fe65+ 
                  ccode_fe66+ ccode_fe67+ ccode_fe68+ ccode_fe69+ ccode_fe70+ ccode_fe71+ ccode_fe72+ ccode_fe73+ 
                  ccode_fe74+ ccode_fe75+ ccode_fe76+ ccode_fe77+ ccode_fe78+ ccode_fe79+ ccode_fe80+ ccode_fe81+ 
                  ccode_fe82+ ccode_fe83+ ccode_fe85+ ccode_fe86+ ccode_fe87+ ccode_fe88+ ccode_fe92+ ccode_fe93+ 
                  ccode_fe94+ ccode_fe95+ ccode_fe96+ ccode_fe97+ 
                  ccode_fe98+ ccode_fe99+ ccode_fe100+ ccode_fe101+ ccode_fe102+ ccode_fe103+ ccode_fe104+ ccode_fe105+ 
                  ccode_fe106+ ccode_fe107+ ccode_fe108+ ccode_fe109+ ccode_fe110+ ccode_fe111+ ccode_fe112+ ccode_fe113+ 
                  ccode_fe114+ ccode_fe115+ ccode_fe116+ ccode_fe117+ ccode_fe118+ ccode_fe119+ ccode_fe120+ ccode_fe121+ 
                  ccode_fe122+ ccode_fe123+ ccode_fe124+ ccode_fe125+ ccode_fe126+ ccode_fe127+ ccode_fe128+ ccode_fe129+ 
                  ccode_fe130+ ccode_fe131+ ccode_fe132+ ccode_fe133+ ccode_fe134+ ccode_fe135+ ccode_fe136+ ccode_fe137+ 
                  ccode_fe138+ y_fe1+ y_fe2+ y_fe3+ y_fe5+ y_fe6+ 
                  y_fe7+ y_fe8+ y_fe9+ y_fe10+ y_fe11+ y_fe12+ y_fe13+ y_fe14+ y_fe15+ y_fe16+ y_fe17+ y_fe18+ y_fe19+ y_fe20+ 
                  y_fe21+ y_fe22+ y_fe23+ y_fe24+ y_fe25+ y_fe26+ y_fe27+ y_fe28+ y_fe29+ y_fe30+ y_fe31+ y_fe32+ y_fe33+ y_fe34+ 
                  y_fe35+ y_fe36+ y_fe37+ y_fe38+ y_fe39+ y_fe40+ y_fe41+ y_fe42+ y_fe43+ y_fe44+ y_fe45+ y_fe46+ y_fe47+ y_fe48+ 
                  y_fe49+ y_fe50+ y_fe51+ y_fe52+ y_fe53+ y_fe54+ y_fe55+ y_fe56+ y_fe57+ y_fe58+ y_fe59+ y_fe60+ y_fe61"


model_sar <- milm(fml=mSTAR, midata=newdata)

model5<- data.frame(model_sar$terms, model_sar$beta, model_sar$SE, model_sar$beta/model_sar$SE)
colnames(model5) <- c("Variable","Est. Coef.","SE","z-score")
print(paste("Model Output"))
print(model5)
